function [ y,b1,mask,Weighting,SG,time_average] = GROG_reconstruction_Preprocessing_N_ref_demo(kdata_corrected_ref,rays_per_frame,nt1)
    [nr,np,nt,nc]=size(kdata_corrected_ref(:,:,:,:));
    kdata_ref = reshape(kdata_corrected_ref(:,:,:,:),[nr,np*nt,nc]);
    kdata = kdata_ref;
    kdata1=kdata;
    [nx,ntviews,nz,nc]=size(kdata1);
    bas=size(kdata_corrected_ref,1);% the 2x full FOV is 512 and we only reconstruct 360 here.
    
    % for zz=1:nz % loop through all slices
    zz=1; % For demonstrate purpose, reconstruct 1 slice only
    kdata2(:,:,1,:) = kdata1;
    kdata=kdata2(:,:,zz,:);
    [nx,ntviews,~,nc]=size(kdata);
    clear kdata2; clear kdata_ref; clear kdata_corrected_ref; clear kdata1 

    %Generating trajectory (For GROG only, NOT for general gridding)
    Traj=Trajectory_GoldenAngle_GROG_N(ntviews,nx,1);%trajectory
    Traj=Traj*(bas/nx);
    %GROG weights calibration
    [Gx,Gy] = GROG.get_Gx_Gy(kdata,Traj);
    
    %Coil sensitivities estimation
    G = GROG.init(kdata,Traj,Gx,Gy,0); % GROG initialization 
    kref = GROG.interp(kdata,G,1,0); % GROG interpolation 
    kref=CropImg(kref,bas,bas);
    ref=gather(squeeze(ifft2c_mri(kref)));
    [b1,time_average]=adapt_array_2d_GROGGRASP(ref);
    b1=single(b1/max(abs(b1(:))));
    
    nline=rays_per_frame;
    nt = floor(ntviews/nline);
    Nqu=bas*pi/2; % Nyquist spokes, please refer to the GROG-GRASP paper for more details about this.
   
    if size(kdata,2) >= ceil(Nqu) 
        G = GROG.init(kdata(:,1:ceil(Nqu),:,:),Traj(:,1:ceil(Nqu)),Gx,Gy,0);
    else
        G = GROG.init(kdata(:,:,:,:),Traj(:,:),Gx,Gy,0);
    end
    D=reshape(G.weight,[sqrt(size(G.weight,1)),sqrt(size(G.weight,1))]);
    DCF=CropImg(D,bas,bas);
    DCF=repmat(DCF,[1,1,nt,nc]);
    clear kdata1; clear kref; clear D;
    
    %% Data sorting
    nspokes=nline;
    clear kdatau_ref,clear ku, clear wu,clear kdatau_venc
    clear Traj_ref;
    for ii=1:nt
        kdatau_ref(:,:,:,ii)=kdata(:,(ii-1)*nspokes+1:ii*nspokes,:);
        Traj_ref (:,:,ii)= Traj(:,(ii-1)*nspokes+1:ii*nspokes);
    end
    clear kdata; clear Traj;
%%
    [nx,ntviews,nc,nt] = size(kdatau_ref);
    kdatau_ref2 = permute(kdatau_ref,[1,2,4,3]);
    clear DCF_U; clear kdatau_ref;
    
    parfor ii = 1:nt1
        G = GROG.init(kdatau_ref2(:,:,ii,:),Traj_ref(:,:,ii),Gx,Gy,0);
        DCF_U(:,:,ii)=reshape(G.weight,[sqrt(size(G.weight,1)),sqrt(size(G.weight,1))]);
        kdata3(:,:,ii,:) = GROG.interp(kdatau_ref2(:,:,ii,:),G,1,0);
    end
    clear kdatau_ref2; clear Traj_ref;

    DCF_U2=repmat(DCF_U,[1,1,1,nc]);
    clear DCF_U; DCF_U = DCF_U2; clear DCF_U2;
    DCF_U= CropImg(DCF_U,bas,bas);
    DCF_U(find(DCF_U==0))=1;
    DCF = DCF(:,:,1:nt1,:);
    Weighting=single(DCF)./single(DCF_U);
    Weighting = Weighting;
    clear DCF; clear DCF_U;
    SoftWeight=single(zeros(bas,bas,nt1,nc));
    SoftWeight(:,:,:,:,1)=(1/4)^(1-1);
    mask=single(kdata3~=0);
    SG=SoftWeight;
    y=kdata3.*sqrt(Weighting);
